##################################################################################################
######################################### READ ME ################################################
##################################################################################################

## Please make sure to save JPART_replication_Markovmodel_JAGS.txt and JPART_replication_data.csv
## in your working directory.

## The replication of the article's Markov model may take several hours depending on your machine.

## Despite the set.seed() command, an element of randomness remains (different machines, different
## package versions, etc.) when using advanced machine learning techniques on relatively small
## data sets. As a result, marginal (less than a decimal) differences may occur between the
## coefficients produced by your replication and those reported in the article.

## Variable overview:
## ID      = Identifier (MS + Year)
## MS      = Member state (country)
## Year    = t, in half-years
## EDP     = Indicates whether MS is under EDP at t
## Vote    = Voting power of MS in Council of the EU
## Govdebt = Gross MS government debt, as percentage of GDP
## Govdef  = Cyclically adjusted MS government deficit, as percentage of GDP
## Govlr   = MS government position on economic left-right dimension
## Goveu   = Ms government position on EU dimension
## Nxtele  = t (half-years) until next parliamentary election in MS
## FN      = Share of MS citizens holding fairly negative opinion of EU
## VN      = Share of MS citizens holding very negative opinion of EU
## ES_MS   = Public Euroscepticism in MS (created from FN + VN below)
## ES_Par  = Euroscepticism in MS parliament
## ES_Cred = Public Euroscepticism in creditor MS (created from ES_MS below)

## All other code (e.g. to replicate PPCs and figures) available from corresponding author on request.

##################################################################################################
######################################### Markov Model ###########################################
##################################################################################################
setwd("YOUR WORKING DIRECTORY")
getwd()

## Loading packages
library(plyr)
library(ggplot2)
library(GGally)
library(lme4)
library(lmerTest)
library(sjstats)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(lattice)
library(arm)
library(brms)
library(optimx)
library(boot)
library(tidyverse)
library(rjags)
library(coda)
library(bayesplot)

## Loading in data
EDP <- read.csv2("JPART_replication_data.csv", header = T, stringsAsFactors = F, na.strings = c('', 'NA'))

## Creating single Euroscepticism variable from VN and FN
EDP$ES_MS <- EDP$VN + EDP$FN

## Creating Euroscepticism variable for creditor states
EDP$ES_Cred <- ((EDP$ES_MS[EDP$MS == 'FIN'] + EDP$ES_MS[EDP$MS == 'NL'] + EDP$ES_MS[EDP$MS == 'GER'] +
                      EDP$ES_MS[EDP$MS == 'AUS'])/4)
EDP$ES_Cred[EDP$MS == 'FIN'] <- ((EDP$ES_MS[EDP$MS == 'NL'] + EDP$ES_MS[EDP$MS == 'GER'] +
                                            EDP$ES_MS[EDP$MS == 'AUS'])/3)
EDP$ES_Cred[EDP$MS == 'NL'] <- ((EDP$ES_MS[EDP$MS == 'FIN'] + EDP$ES_MS[EDP$MS == 'GER'] +
                                           EDP$ES_MS[EDP$MS == 'AUS'])/3)
EDP$ES_Cred[EDP$MS == 'GER'] <- ((EDP$ES_MS[EDP$MS == 'NL'] + EDP$ES_MS[EDP$MS == 'FIN'] +
                                            EDP$ES_MS[EDP$MS == 'AUS'])/3)
EDP$ES_Cred[EDP$MS == 'AUS'] <- ((EDP$ES_MS[EDP$MS == 'FIN'] + EDP$ES_MS[EDP$MS == 'GER'] +
                                            EDP$ES_MS[EDP$MS == 'NL'])/3)

## MS to numeric
EDP$MS <- factor(EDP$MS)
EDP$MS <- droplevels(EDP$MS)
EDP$MS <- as.integer(EDP$MS)

## Creating required functions:
twoSD <- function (x) {
  (x - mean(x, na.rm = TRUE))/(2 * sd(x, na.rm = TRUE))}   ## Standardizing by 2sd & grand mean centering

getSE <- function (x) {                                    ## shows SE
  sqrt(diag(vcov(x)))}

meanc <- function (x) {                                    ## Normal grand mean centering
  (x - mean(x, na.rm = TRUE))}

as.numeric.factor <- function(x) {as.numeric(levels(x))[x]}

## Group mean centering predictor variables
EDP$Vote_gc             <- meanc(EDP$Vote)
EDP$ES_MS_gc             <- meanc(EDP$ES_MS)
EDP$GovDebt_gc           <- meanc(EDP$GovDebt)
EDP$GovDef_gc            <- meanc(EDP$GovDef)
EDP$Govlr_gc            <- meanc(EDP$Govlr)
EDP$Goveu_gc            <- meanc(EDP$Goveu)
EDP$ES_Par_gc            <- meanc(EDP$ES_Par)
EDP$Nxtele_gc         <- meanc(EDP$Nxtele)
EDP$ES_Cred_gc         <- meanc(EDP$ES_Cred)

## and standardizing non-binary by two SD
EDP$Vote_gc             <- EDP$Vote_gc/(2*sd(EDP$Vote, na.rm = T))
EDP$ES_MS_gc             <- EDP$ES_MS_gc/(2*sd(EDP$ES_MS, na.rm = T))
EDP$GovDebt_gc           <- EDP$GovDebt_gc/(2*sd(EDP$GovDebt, na.rm = T))
EDP$GovDef_gc            <- EDP$GovDef_gc/(2*sd(EDP$GovDef, na.rm = T))
EDP$Govlr_gc            <- EDP$Govlr_gc/(2*sd(EDP$Govlr, na.rm = T))
EDP$Goveu_gc            <- EDP$Goveu_gc/(2*sd(EDP$Goveu, na.rm = T))
EDP$ES_Par_gc            <- EDP$ES_Par_gc/(2*sd(EDP$ES_Par, na.rm = T))
EDP$Nxtele_gc         <- EDP$Nxtele_gc/(2*sd(EDP$Nxtele, na.rm = T))
EDP$ES_Cred_gc         <- EDP$ES_Cred_gc/(2*sd(EDP$ES_Cred, na.rm = T))

## Formatting structure: JAGS requires all variables to be entered as seperate matrices
keep <- c("MS", "Year", "EDP")    #"GovDebt_gc", "GovDef_gc", "EMU_gc", "Vote_gc","Goveu_gc", "Govlr_gc", "ES_Par_gc", "Nxtele_gc", "ES_MS_gc")
state <- EDP[,keep]
state <- reshape(state, idvar = "MS", timevar = "Year", direction = "wide")

colnames(state) <- c("MS", c(unique(EDP$Year)))
state <- as.data.frame(t(state))
colnames(state) <- as.character(unlist(state[1,]))
state <- state[-1, ]
data <- state
state <- state+1

ES_MS <- EDP[, c("MS", "Year", "ES_MS_gc")]
ES_MS <- reshape(ES_MS, idvar = "MS", timevar = "Year", direction = "wide")
colnames(ES_MS) <- c("MS", c(unique(EDP$Year)))
ES_MS <- as.data.frame(t(ES_MS))
colnames(ES_MS) <- as.character(unlist(ES_MS[1,]))
ES_MS <- ES_MS[-1, ]

GovDebt <- EDP[, c("MS", "Year", "GovDebt_gc")]
GovDebt <- reshape(GovDebt, idvar = "MS", timevar = "Year", direction = "wide")
colnames(GovDebt) <- c("MS", c(unique(EDP$Year)))
GovDebt <- as.data.frame(t(GovDebt))
colnames(GovDebt) <- as.character(unlist(GovDebt[1,]))
GovDebt <- GovDebt[-1, ]

GovDef <- EDP[, c("MS", "Year", "GovDef_gc")]
GovDef <- reshape(GovDef, idvar = "MS", timevar = "Year", direction = "wide")
colnames(GovDef) <- c("MS", c(unique(EDP$Year)))
GovDef <- as.data.frame(t(GovDef))
colnames(GovDef) <- as.character(unlist(GovDef[1,]))
GovDef <- GovDef[-1, ]

Vote <- EDP[, c("MS", "Year", "Vote_gc")]
Vote <- reshape(Vote, idvar = "MS", timevar = "Year", direction = "wide")
colnames(Vote) <- c("MS", c(unique(EDP$Year)))
Vote <- as.data.frame(t(Vote))
colnames(Vote) <- as.character(unlist(Vote[1,]))
Vote <- Vote[-1, ]

Govlr <- EDP[, c("MS", "Year", "Govlr_gc")]
Govlr <- reshape(Govlr, idvar = "MS", timevar = "Year", direction = "wide")
colnames(Govlr) <- c("MS", c(unique(EDP$Year)))
Govlr <- as.data.frame(t(Govlr))
colnames(Govlr) <- as.character(unlist(Govlr[1,]))
Govlr <- Govlr[-1, ]

Goveu <- EDP[, c("MS", "Year", "Goveu_gc")]
Goveu <- reshape(Goveu, idvar = "MS", timevar = "Year", direction = "wide")
colnames(Goveu) <- c("MS", c(unique(EDP$Year)))
Goveu <- as.data.frame(t(Goveu))
colnames(Goveu) <- as.character(unlist(Goveu[1,]))
Goveu <- Goveu[-1, ]

ES_Par <- EDP[, c("MS", "Year", "ES_Par_gc")]
ES_Par <- reshape(ES_Par, idvar = "MS", timevar = "Year", direction = "wide")
colnames(ES_Par) <- c("MS", c(unique(EDP$Year)))
ES_Par <- as.data.frame(t(ES_Par))
colnames(ES_Par) <- as.character(unlist(ES_Par[1,]))
ES_Par <- ES_Par[-1, ]

Nxtele <- EDP[, c("MS", "Year", "Nxtele_gc")]
Nxtele <- reshape(Nxtele, idvar = "MS", timevar = "Year", direction = "wide")
colnames(Nxtele) <- c("MS", c(unique(EDP$Year)))
Nxtele <- as.data.frame(t(Nxtele))
colnames(Nxtele) <- as.character(unlist(Nxtele[1,]))
Nxtele <- Nxtele[-1, ]

ES_Cred <- EDP[, c("MS", "Year", "ES_Cred_gc")]
ES_Cred <- reshape(ES_Cred, idvar = "MS", timevar = "Year", direction = "wide")
colnames(ES_Cred) <- c("MS", c(unique(EDP$Year)))
ES_Cred <- as.data.frame(t(ES_Cred))
colnames(ES_Cred) <- as.character(unlist(ES_Cred[1,]))
ES_Cred <- ES_Cred[-1, ]

# Specify initial values if needed
beta.inits <- matrix(data= c(NA, -0.3,0.3, NA),  # no inits for the diagonal, because the diagonals are not estimated, but fixed at 0.
                     nrow=2, ncol=2, byrow=T)

# Use jags.model function to create a 'jags' object
require(rjags)
K <- 2
M <- 2
T <- 27
N <- 28
print(date())

## JAGS model
model <- jags.model(file="JPART_replication_Markovmodel_JAGS.txt",
                    data=list(m=state, GovDebt=GovDebt, GovDef=GovDef, Vote=Vote, ES_Cred=ES_Cred, ES_MS=ES_MS, Goveu=Goveu,
                              ES_Par=ES_Par, Govlr=Govlr, Nxtele=Nxtele,
                              statealpha=rep(1,times=M), M=M, T=T, N=N),
                    #    inits=list(beta1=beta.inits,beta2=beta.inits,beta3=beta.inits,beta11=beta.inits,beta5=beta.inits,
                    #               beta6=beta.inits,beta7=beta.inits,beta8=beta.inits,beta9=beta.inits,beta10=beta.inits,beta12=beta.inits),
                    n.chains=4, n.adapt=1000, quiet=F)

print(date())

## Adaptation of the sampling algorithms
set.seed(113)
repeat{ adapted <- adapt(model, n.iter=1000, end.adaptation=FALSE)    ; 
if(adapted == TRUE)  { break } } # if adaptation is finished the loop ends
adapt(model, n.iter=10, end.adaptation=TRUE)

## For a quicker yet dirtier replication, reduce the number of iterations (burnin and stored),
## or the number of chains (n.chains, line 193)
burnin <- (10 * 1000)
stored <- (50 * 1000)
thinparameter <- 10

## Do burnin iterations, then set monitors and collect samples
print(date())  ;  update(model, n.iter=burnin)  ;  print(date())
opslaan <- c('alpha', 'delta', 'beta1', 'beta2', 'beta3', 'beta4', 'beta5', 'beta6', 'beta7', 'beta8',
             'beta9', 'beta10', 'beta11', 'beta12',
             'm_pred', 'alpha.means', 'alpha.taus', 'alpha.sigmasq12', 'alpha.Sigma', 'alpha.sigmasq21',
             'means', 'sigma2e', 'sigma_e', 'probs', 'u', 'umean', 'usigma', 'utau', 
             'usigmasq', 'taus', 'probs1')
resultaten <- coda.samples(model, variable.names=opslaan, n.iter=stored, 
                           thin=thinparameter)  ;  print(date())

## Store summaries of the posterior distributions.
## All OMM coefficients can be extracted from this object.
estimates <- summary(resultaten, quantiles = c(0.025, 0.5, 0.975))

## Highest Posterior Density Intervals
require(coda) 
combinedresults <- as.mcmc(do.call(rbind,resultaten)) 
HPDs <- HPDinterval(combinedresults, prob = 0.95)
HPDs <- round(HPDs, digits = 4)

## Posterior estimates for verification
posterior12 <- mcmc_intervals(resultaten, pars = rev(c("alpha.means[1]", "beta1[1,2]", "beta2[1,2]", "beta3[1,2]", "beta4[1,2]",
                                                       "beta5[1,2]", "beta6[1,2]", "beta7[1,2]", "beta8[1,2]",
                                                       "beta9[1,2]","beta10[1,2]", "beta11[1,2]", "beta12[1,2]")),
                              point_est = 'median') + scale_y_discrete(labels = rev(c('Intercept', 'Gross Gov. Debt', 'Cycl. Adj. Deficit',
                                                                                      'Voting power', 'Euroscepticism Creditor MS', 'Public Euroscepticism MS', 'Mobilization',
                                                                                      'Gov. Pos. EU', 'Gov. Pos. Left-Right', 'Electoral cycle',
                                                                                      'Eurosc. MS x Deficit', 'Eurosc. MS x Mobilization', 
                                                                                      'Eurosc. MS x Eurosc. Cred'))) +
  ggtitle('OMM results for entering EDP') + theme_classic()
posterior12

posterior21 <- mcmc_intervals(resultaten, pars = rev(c("alpha.means[2]", "beta1[2,1]", "beta2[2,1]", "beta3[2,1]", "beta4[2,1]",
                                                       "beta5[2,1]", "beta6[2,1]", "beta7[2,1]", "beta8[2,1]",
                                                       "beta9[2,1]","beta10[2,1]", "beta11[2,1]", "beta12[2,1]")),
                              point_est = 'median') + scale_y_discrete(labels = rev(c('Intercept', 'Gross Gov. Debt', 'Cycl. Adj. Deficit',
                                                                                      'Voting power', 'Euroscepticism Creditor MS', 'Public Euroscepticism MS', 'Mobilization',
                                                                                      'Gov. Pos. EU', 'Gov. Pos. Left-Right', 'Electoral cycle',
                                                                                      'Eurosc. MS x Deficit', 'Eurosc. MS x Mobilization', 
                                                                                      'Eurosc. MS x Eurosc. Cred')))  +
  ggtitle('OMM results for leaving EDP') + theme_classic()
posterior21


#################################################################################################################
################################### Simpler models (non-bayes) ##################################################
#################################################################################################################
EDP2 <- EDP
colnames(EDP2)[4] <- "state"
EDP2 <- EDP2 %>%
  group_by(MS) %>%
  mutate(enter = as.numeric(state != dplyr::lag(state)))
EDP$enter <- EDP2$enter
rm(EDP2)
EDP$enter[is.na(EDP$enter)] <- 0
EDP$enter[EDP$enter == 1 & EDP$EDP == 0] <- 0
EDP3 <- EDP[,c('enter', 'EDP', 'ID')]
EDP2 <- EDP[!(EDP$enter == 0 & EDP$EDP == 1),]

## Models 0-4
model0 <- glmer(EDP ~ 1 + (1|MS), data = EDP2,
                family=binomial(link="logit"))

model4 <- glmer(EDP ~ 1 + GovDebt_gc + GovDef_gc + Vote_gc + ES_Cred_gc + ES_MS_gc +
                  ES_Par_gc + Goveu_gc + Govlr_gc + Nxtele_gc +
                  ES_MS_gc*GovDef_gc + ES_MS_gc*ES_Par_gc + ES_MS_gc*ES_Cred_gc +
                  (1 | MS) + (1|Year), data = EDP2,
                family=binomial(link="logit"), glmerControl(optimizer = 'bobyqa'))

model3 <- glmer(EDP ~ 1 + GovDebt_gc + GovDef_gc + Vote_gc + ES_Cred_gc + ES_MS_gc +
                  ES_Par_gc + Goveu_gc + Govlr_gc + Nxtele_gc +
                  ES_MS_gc*GovDef_gc + ES_MS_gc*ES_Par_gc + ES_MS_gc*ES_Cred_gc +
                  (1 | MS), data = EDP2,
                family=binomial(link="logit"), glmerControl(optimizer = 'bobyqa'))

model2 <- glmer(EDP ~ 1 + GovDebt_gc + GovDef_gc + Vote_gc + Goveu_gc + Govlr_gc + Nxtele_gc +
                  (1 | MS), data = EDP2,
                family=binomial(link="logit"), glmerControl(optimizer = 'bobyqa'))

model1 <- glmer(EDP ~ 1 + GovDebt_gc + GovDef_gc +
                  (1 | MS), data = EDP2,
                family=binomial(link="logit"), glmerControl(optimizer = 'bobyqa'))

## Table 4 as reported in appendix
tab_model(model0, model1, model2, model3, model4, transform = NULL)

## Table 5 as reported in appendix
ano <- anova(model0, model1, model2, model3, model4, test = "Chisq")



#################################################################################################################
################################### Bayesian Survival model #####################################################
#################################################################################################################
EDP2 <- as.data.frame(EDP2 %>%
                        group_by(MS) %>%
                        mutate(enter1 = lag(enter)))
EDP2$enter1[is.na(EDP2$enter1)] <- 0

EDP3 <- EDP2 %>%
  group_by(MS) %>%
  mutate(t = ave(enter1, cumsum(enter1 == 1), FUN = seq))

library('survival')
library('survminer')
library('splines2')
surv <- brm(t | cens(1-enter) ~ 1 + GovDebt_gc + GovDef_gc + Vote_gc + ES_Cred_gc + ES_MS_gc +
              ES_Par_gc + Goveu_gc + Govlr_gc + Nxtele_gc +
              ES_MS_gc*GovDef_gc + ES_MS_gc*ES_Par_gc + ES_MS_gc*ES_Cred_gc +
              (1|MS),
            data = EDP3,
            family  = brmsfamily("cox"),
            warmup  = 2000,
            iter    = 5000,
            control = list(adapt_delta = .99))
model5 <- summary(surv)
model5 <- as.data.frame(model5$fixed)
model5$hazrat <- exp(model5$Estimate)
model5 <- round(model5, 3)
model5

